Laboratorio (hands-on): Reconstrucción básica de una radiografía (Beer–Lambert + haz paralelo)

Fuente → Objeto → Detector → (normalización + log) → Radiografía lista para ver

Autor/a

Ph.D. Pablo Eduardo Caicedo Rodríguez

Fecha de publicación

16 de febrero de 2026

1) Introducción

En este laboratorio vas a construir una “radiografía” sintética desde cero, entendiendo la cadena fuente → objeto → detector. La idea central es simple: si un rayo atraviesa más material (o material más atenuante), llegan menos fotones al detector.

¿Qué entra y qué sale?

Entra (fuente): cuántos fotones salen (lo llamaremos I0).
Entra (objeto): un mapa 2D mu[y,x] que dice qué tan atenuante es cada zona.
Sale (detector): cuántos fotones llegan por rayo (lo llamaremos I_det).
Sale (radiografía): una imagen final rad = -log(I_det/I0) que es proporcional a “cuánta atenuación acumuló” cada rayo.

Objetivos de aprendizaje (al finalizar deberías poder)

Explicar con palabras simples Beer–Lambert.
Construir phantoms 2D con regiones de distintos materiales.
Simular proyecciones con haz paralelo (sumas en filas).
Interpretar I_det (más fotones vs menos fotones) y rad (más atenuación vs menos).

Mapa del pipeline (una sola vista)

I0 (fuente)  →  tau (suma de μ en el camino)  →  I_det = I0·exp(-tau)  →  rad = -log(I_det/I0)

Glosario rápido (lo mínimo para no perderse)

Símbolo/variable	Qué es	Cómo se interpreta
`mu` (μ)	Mapa 2D del objeto	μ grande = material más atenuante
`tau` (τ)	Atenuación acumulada por rayo	τ grande = el rayo vio “más bloqueo”
`I0`	Intensidad sin objeto	referencia (brillo de la fuente)
`I_det` (I)	Intensidad medida con objeto	baja cuando τ sube
`rad`	`-log(I_det/I0)`	sube cuando τ sube (radiografía)

Nota

Regla mental que usarás todo el tiempo Más mu (o más “recorrido” dentro del objeto) → mayor atenuación acumulada tau → menor intensidad en el detector I_det → mayor valor en la radiografía rad (después del log).

2) Preparación (imports y parámetros)

Trabajaremos con arreglos 2D N×N. En 2D, el detector “natural” es una señal 1D (un valor por fila). Para visualizarlo como imagen, repetiremos esa señal (tiling) y la mostraremos como 2D.

Código

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["figure.dpi"] = 120

N = 256       # tamaño del phantom (rápido en laptop)
dx = 1.0      # tamaño de pixel (unidades arbitrarias)
eps = 1e-12   # evita log(0) si hay valores muy pequeños

3) Modelo físico: Beer–Lambert (muy claro y sin extras)

3.1 La idea (en palabras)

Imagina el objeto como una cortina entre la fuente y el detector. - Cortina “ligera” (μ pequeño) → pasan más fotones. - Cortina “pesada” (μ grande) → pasan menos fotones.

μ(x,y) es un número que indica qué tan “bloqueador” es el material en ese punto. Un rayo no ve un solo punto: atraviesa una ruta y va acumulando bloqueo.

3.2 La ecuación

\[ I = I_0 \, \exp\left(-\int \mu(x,y)\,dl\right) \]

Traducción directa (sin misterio): - I0: fotones que llegarían si NO hubiera objeto. - I: fotones que llegan DESPUÉS de atravesar el objeto. - ∫ μ dl (lo llamaremos tau): “bloqueo total” acumulado a lo largo del camino del rayo.

3.3 Por qué luego usamos log

Si definimos

\[ \tau = \int \mu\,dl, \]

entonces

\[ \frac{I}{I_0} = e^{-\tau} \quad\Rightarrow\quad -\log\left(\frac{I}{I_0}\right)=\tau. \]

O sea: el log “deshace” la exponencial y nos devuelve algo proporcional a la atenuación acumulada. Eso es lo que aquí llamamos “radiografía reconstruida”.

Advertencia

Qué NO incluye este laboratorio (a propósito): no hay scatter, no hay ruido Poisson, no hay endurecimiento de haz, no hay blur del detector. Esto es para aprender el esqueleto del pipeline sin distracciones.

4) Construcción de un phantom (mapa 2D de materiales μ)

4.1 Qué es un phantom (didáctico)

Un phantom es un “objeto inventado” para practicar. Aquí es una imagen 2D de números mu[y,x]: cada pixel dice “qué material hay” (qué tanto atenúa). - mu pequeño → zona más “transparente” (deja pasar más fotones). - mu grande → zona más “opaca” (deja pasar menos fotones).

El fondo representa un material base (mu_fondo). Las formas internas tienen otro μ para crear contraste.

4.2 Cómo “dibujamos” materiales con máscaras

En NumPy podemos construir una máscara booleana (True/False). - True = “píntalo con el material del círculo/bloque”. - False = “déjalo como fondo”.

En el siguiente bloque creamos dos phantoms: A) fondo + círculo, B) fondo + bloque + círculo.

Código

def crear_phantom_circulo(N, mu_fondo=0.02, mu_circulo=0.06, radio_frac=0.25, centro=None):
    """Phantom A: fondo uniforme + un círculo."""
    if centro is None:
        centro = (N//2, N//2)  # (y, x)

    # Grid de coordenadas
    y = np.arange(N)[:, None]
    x = np.arange(N)[None, :]

    cy, cx = centro
    r = radio_frac * N

    # 1) Empieza con fondo
    phantom = np.full((N, N), mu_fondo, dtype=float)

    # 2) Máscara del círculo (True dentro del círculo)
    mascara = (x - cx)**2 + (y - cy)**2 <= r**2

    # 3) "Pinta" el círculo cambiando mu
    phantom[mascara] = mu_circulo
    return phantom

def crear_phantom_dos_materiales(N, mu_fondo=0.02, mu_bloque=0.10, mu_circulo=0.05):
    """Phantom B: fondo + bloque (más atenuante) + círculo (intermedio)."""
    phantom = np.full((N, N), mu_fondo, dtype=float)

    # Bloque rectangular (material más atenuante)
    y1, y2 = int(0.35*N), int(0.75*N)
    x1, x2 = int(0.15*N), int(0.45*N)
    phantom[y1:y2, x1:x2] = mu_bloque

    # Círculo (material intermedio)
    y = np.arange(N)[:, None]
    x = np.arange(N)[None, :]
    cy, cx = int(0.45*N), int(0.70*N)
    r = int(0.15*N)
    mascara = (x - cx)**2 + (y - cy)**2 <= r**2
    phantom[mascara] = mu_circulo

    return phantom

mu_A = crear_phantom_circulo(N)
mu_B = crear_phantom_dos_materiales(N)

fig, ax = plt.subplots(1, 2, figsize=(9, 4))

im0 = ax[0].imshow(mu_A, cmap="gray")
ax[0].set_title("Phantom A: círculo en fondo (mapa μ)")
plt.colorbar(im0, ax=ax[0], fraction=0.046)

im1 = ax[1].imshow(mu_B, cmap="gray")
ax[1].set_title("Phantom B: dos materiales (mapa μ)")
plt.colorbar(im1, ax=ax[1], fraction=0.046)

for a in ax:
    a.set_axis_off()

plt.tight_layout()
plt.show()

Tip

Lectura rápida del phantom En la imagen del phantom NO estás viendo una radiografía. Estás viendo “materiales”: dónde μ es mayor, ese material atenúa más.

Ejercicio guiado 1 (predicción antes de correr el pipeline)

Si aumentas mu_circulo, ¿la intensidad I_det en esas filas sube o baja? ¿por qué?
Si aumentas mu_bloque, ¿qué esperas que pase con I_det y con rad?

Escribe tu predicción en 2–3 líneas. Luego la comprobarás cuando veas I_det y rad.

5) Generador polienergético idealizado (espectro simple)

5.1 Qué significa “polienergético” aquí

En rayos X reales, la fuente emite fotones con diferentes energías (un espectro). Aquí lo simplificamos a 3 energías discretas con pesos. Esto es didáctico.

Importante: en este laboratorio NO hacemos que μ dependa de energía. Por eso, el “espectro” se reduce a un I0_total efectivo, pero nos sirve para entender la idea de “sumar contribuciones”.

5.2 Implementación (3 energías + pesos)

Código

energias_keV = np.array([30, 60, 90], dtype=float)
pesos = np.array([0.50, 0.30, 0.20], dtype=float)
pesos = pesos / pesos.sum()

I0_total = 1e5
I0_por_energia = I0_total * pesos

print("Energías (keV):", energias_keV)
print("Pesos:", pesos)
print("I0 por energía:", I0_por_energia)
print("I0 efectivo:", I0_por_energia.sum())

# Visualización del espectro discreto (solo para ver pesos)
fig, ax = plt.subplots(figsize=(5.5, 3.2))
ax.bar(energias_keV.astype(int), pesos)
ax.set_title("Espectro simple (discreto): pesos por energía")
ax.set_xlabel("Energía (keV)")
ax.set_ylabel("Peso (suma = 1)")
plt.tight_layout()
plt.show()

Energías (keV): [30. 60. 90.]
Pesos: [0.5 0.3 0.2]
I0 por energía: [50000. 30000. 20000.]
I0 efectivo: 100000.0

Pregunta corta (1 línea): Si μ no cambia con energía, ¿por qué esperarías que rad = -log(I/I0) cambie muy poco al cambiar pesos?

6) Simulación de proyecciones (haz paralelo, integral discreta)

6.1 Qué es “proyectar” en este laboratorio

Con haz paralelo y en 2D, usaremos rayos horizontales. Para cada fila y, el rayo atraviesa todos los pixeles x y acumula atenuación:

\[ \tau(y) \approx \sum_x \mu(y,x)\,\Delta x \]

tau(y) es una señal 1D: un valor por fila (un valor por rayo).

6.2 Implementación (sumar en x)

Código

def proyectar_haz_paralelo(mu, dx=1.0):
    """Proyección 2D→1D: tau[y] = sum_x mu[y,x]*dx"""
    tau = np.sum(mu, axis=1) * dx
    return tau

tau_A = proyectar_haz_paralelo(mu_A, dx=dx)
tau_B = proyectar_haz_paralelo(mu_B, dx=dx)

fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].plot(tau_A)
ax[0].set_title("Camino óptico τ(y) - A")
ax[0].set_xlabel("índice y (pixel)")
ax[0].set_ylabel("τ (aprox.)")

ax[1].plot(tau_B)
ax[1].set_title("Camino óptico τ(y) - B")
ax[1].set_xlabel("índice y (pixel)")

plt.tight_layout()
plt.show()

6.2b Conectar “una fila del phantom” con un valor de τ

Para que no se vuelva abstracto: el rayo de la fila y0 atraviesa todos los pixeles de esa fila. Entonces tau[y0] es literalmente la suma de mu[y0, :] (por dx). En el siguiente bloque miramos una fila y verificamos la suma.

Código

# Escoge una fila (rayo) para inspeccionar
y0 = N // 2

# Perfil de mu en esa fila
perfil_A = mu_A[y0, :]
perfil_B = mu_B[y0, :]

print(f"Fila y0={y0}: tau_A[y0]={tau_A[y0]:.4f}  (suma discreta)")
print(f"Fila y0={y0}: tau_B[y0]={tau_B[y0]:.4f}  (suma discreta)")

fig, ax = plt.subplots(1, 2, figsize=(9, 3))
ax[0].plot(perfil_A)
ax[0].set_title("Perfil μ en una fila (A)")
ax[0].set_xlabel("x (pixel)")
ax[0].set_ylabel("μ")

ax[1].plot(perfil_B)
ax[1].set_title("Perfil μ en una fila (B)")
ax[1].set_xlabel("x (pixel)")

plt.tight_layout()
plt.show()

Fila y0=128: tau_A[y0]=10.2800  (suma discreta)
Fila y0=128: tau_B[y0]=13.4100  (suma discreta)

Nota

Interpretación de τ(y) Si en cierta fila el rayo cruza el círculo o el bloque, τ(y) sube (porque suma μ más grande). Más τ → menos intensidad en el detector (por la exponencial).

6.3 Para ver 1D como “imagen” (solo visual)

Esto NO agrega física: solo repite la señal 1D en columnas para verla como imagen.

Código

def tile_a_imagen(signal_1d, N):
    """Convierte una señal 1D (N,) a una imagen (N,N) repitiéndola en x."""
    return np.tile(signal_1d[:, None], (1, N))

tau_A_img = tile_a_imagen(tau_A, N)
tau_B_img = tile_a_imagen(tau_B, N)

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
im0 = ax[0].imshow(tau_A_img, cmap="magma")
ax[0].set_title("τ(y) como imagen (tiling) - A")
plt.colorbar(im0, ax=ax[0], fraction=0.046)

im1 = ax[1].imshow(tau_B_img, cmap="magma")
ax[1].set_title("τ(y) como imagen (tiling) - B")
plt.colorbar(im1, ax=ax[1], fraction=0.046)

for a in ax:
    a.set_axis_off()

plt.tight_layout()
plt.show()

7) Detector ideal (sin blur, sin ruido)

7.1 Qué calcula el detector

El detector ideal “cuenta” fotones que llegan por rayo. Para cada energía k:

\[ I_k(y) = I0_k\,e^{-\tau(y)} \]

y luego sumamos todas las energías:

\[ I(y) = \sum_k I_k(y). \]

7.2 Implementación

Código

def detector_ideal_polienergetico(tau, I0_por_energia):
    """I(y) = sum_k I0_k * exp(-tau(y))"""
    I = np.zeros_like(tau, dtype=float)
    for I0_k in I0_por_energia:
        I += I0_k * np.exp(-tau)
    return I

I_det_A = detector_ideal_polienergetico(tau_A, I0_por_energia)
I_det_B = detector_ideal_polienergetico(tau_B, I0_por_energia)

fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].plot(I_det_A)
ax[0].set_title("Intensidad ideal en detector I(y) - A")
ax[0].set_xlabel("índice y (pixel)")
ax[0].set_ylabel("conteos (ideal)")

ax[1].plot(I_det_B)
ax[1].set_title("Intensidad ideal en detector I(y) - B")
ax[1].set_xlabel("índice y (pixel)")

plt.tight_layout()
plt.show()

7.2b Ver la exponencial con tus ojos (τ vs I)

Si Beer–Lambert es correcto, al aumentar tau la intensidad debe caer de forma exponencial. En este bloque graficamos I_det contra tau para ver esa caída.

Código

fig, ax = plt.subplots(figsize=(5.8, 4))
ax.scatter(tau_A, I_det_A, s=10, alpha=0.6, label="A")
ax.scatter(tau_B, I_det_B, s=10, alpha=0.6, label="B")
ax.set_title("Relación exponencial: I_det vs τ")
ax.set_xlabel("τ (atenuación acumulada)")
ax.set_ylabel("I_det (conteos ideales)")
ax.legend()
plt.tight_layout()
plt.show()

Tip

Cómo “leer” I_det Picos altos = llegaron muchos fotones (poca atenuación). Valles bajos = llegaron pocos fotones (mucha atenuación).

Ejercicio guiado 2 (comprobación de tu predicción)

Vuelve al Ejercicio 1. Ahora que ves I_det, confirma si tu predicción fue correcta y escribe 2 líneas: - “Subí μ en X → I_det en esa zona bajó/subió porque …”.

8) Reconstrucción tipo radiografía: normalización + log

8.1 Por qué hacemos estos dos pasos

Normalización: I_det/I0 elimina el “brillo” de la fuente y deja solo el efecto del objeto.
Log negativo: como I/I0 = exp(-tau), entonces -log(I/I0) = tau.

Esto es clave: la radiografía reconstruida rad es proporcional al “camino óptico” tau.

8.2 Implementación

Código

def reconstruir_radiografia(I_det, I0_total, eps=1e-12):
    """rad(y) = -log( I_det(y) / I0_total )"""
    I_norm = (I_det + eps) / (I0_total + eps)
    rad = -np.log(I_norm)
    return rad, I_norm

rad_A, I_norm_A = reconstruir_radiografia(I_det_A, I0_total, eps=eps)
rad_B, I_norm_B = reconstruir_radiografia(I_det_B, I0_total, eps=eps)

# Verificación visual: rad ~ tau
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].plot(tau_A, label="tau (verdad)")
ax[0].plot(rad_A, "--", label="rad (-log)")
ax[0].set_title("A: tau vs rad")
ax[0].legend()

ax[1].plot(tau_B, label="tau (verdad)")
ax[1].plot(rad_B, "--", label="rad (-log)")
ax[1].set_title("B: tau vs rad")
ax[1].legend()

plt.tight_layout()
plt.show()

Nota

Interpretación de rad Valores altos en rad significan: el rayo vio más material (o material más atenuante). Esto suele ser lo que queremos “ver” en una radiografía.

Ejercicio guiado 3 (invarianza a I0)

Multiplica I0_total por 10 (por ejemplo 1e6) y repite el pipeline. - ¿I_det cambia? (sí) - ¿rad cambia mucho? (debería cambiar muy poco)

9) Visualización: ventaneo (contraste simple)

9.1 Por qué ventaneamos

Una radiografía es una matriz de valores, pero el ojo necesita un rango limitado para ver detalles. Ventaneo = escoger un intervalo [low, high] y mapearlo a [0,1]. - Ancho pequeño → más contraste, pero se saturan zonas. - Ancho grande → menos contraste, más tonos visibles.

9.2 Implementación

Código

def ventaneo_lineal(img, centro, ancho):
    """
    Ventaneo simple:
    - recorta a [centro - ancho/2, centro + ancho/2]
    - reescala a [0, 1]
    """
    low = centro - ancho/2
    high = centro + ancho/2
    out = np.clip(img, low, high)
    out = (out - low) / (high - low + 1e-12)
    return out

# Convertimos señales 1D a imágenes 2D (tiling) para visualizar
rad_A_img = tile_a_imagen(rad_A, N)
rad_B_img = tile_a_imagen(rad_B, N)

# Ventaneo automático simple con percentiles
centro_A = float(np.median(rad_A))
ancho_A = float(np.percentile(rad_A, 95) - np.percentile(rad_A, 5))
centro_B = float(np.median(rad_B))
ancho_B = float(np.percentile(rad_B, 95) - np.percentile(rad_B, 5))

rad_A_win = ventaneo_lineal(rad_A_img, centro_A, ancho_A)
rad_B_win = ventaneo_lineal(rad_B_img, centro_B, ancho_B)

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
ax[0].imshow(rad_A_win, cmap="gray", vmin=0, vmax=1)
ax[0].set_title("Radiografía A (ventaneada)")
ax[1].imshow(rad_B_win, cmap="gray", vmin=0, vmax=1)
ax[1].set_title("Radiografía B (ventaneada)")

for a in ax:
    a.set_axis_off()

plt.tight_layout()
plt.show()

9.3 ¿Qué está haciendo el ventaneo? (histograma)

El histograma te muestra cuántos pixeles tienen cierto valor de rad. Al ventanear, elegimos un rango [low, high] y todo lo demás se recorta. Aquí dibujamos ese rango sobre el histograma para entenderlo.

Código

low_A = centro_A - ancho_A/2
high_A = centro_A + ancho_A/2

fig, ax = plt.subplots(figsize=(6.2, 3.6))
ax.hist(rad_A, bins=40)
ax.axvline(low_A, linestyle='--', label='low (A)')
ax.axvline(high_A, linestyle='--', label='high (A)')
ax.set_title("Histograma de rad (A) y rango de ventaneo")
ax.set_xlabel("rad")
ax.set_ylabel("número de filas (pixeles en 1D)")
ax.legend()
plt.tight_layout()
plt.show()

Ejercicio guiado 4 (control de contraste)

Prueba estos cambios y describe en 2–3 líneas qué notas: - Reduce ancho_A a la mitad. ¿Qué detalle aparece y qué se satura? - Aumenta ancho_A al doble. ¿Qué se ve más “plano” y por qué?

10) Comparación final: phantom vs detector vs radiografía

Aquí juntamos todo para leer el pipeline de un vistazo: - mu = “qué hay” (materiales). - I_det = “qué midió el detector” (conteos). - rad = “radiografía” (atenuación acumulada).

Código

def resumen_pipeline(mu, I0_total, I0_por_energia, dx=1.0, eps=1e-12):
    tau = np.sum(mu, axis=1) * dx
    I_det = detector_ideal_polienergetico(tau, I0_por_energia)
    rad, _ = reconstruir_radiografia(I_det, I0_total, eps=eps)
    return tau, I_det, rad

tauA, IA, radA = resumen_pipeline(mu_A, I0_total, I0_por_energia, dx=dx, eps=eps)
tauB, IB, radB = resumen_pipeline(mu_B, I0_total, I0_por_energia, dx=dx, eps=eps)

fig, ax = plt.subplots(2, 3, figsize=(10, 6))

ax[0,0].imshow(mu_A, cmap="gray"); ax[0,0].set_title("A: phantom μ(y,x)")
ax[1,0].imshow(mu_B, cmap="gray"); ax[1,0].set_title("B: phantom μ(y,x)")

ax[0,1].imshow(tile_a_imagen(IA, N), cmap="gray"); ax[0,1].set_title("A: detector I(y) (tiling)")
ax[1,1].imshow(tile_a_imagen(IB, N), cmap="gray"); ax[1,1].set_title("B: detector I(y) (tiling)")

ax[0,2].imshow(tile_a_imagen(radA, N), cmap="inferno"); ax[0,2].set_title("A: rad = -log(I/I0)")
ax[1,2].imshow(tile_a_imagen(radB, N), cmap="inferno"); ax[1,2].set_title("B: rad = -log(I/I0)")

for a in ax.ravel():
    a.set_axis_off()

plt.tight_layout()
plt.show()

11) Mini-proyecto final (entregable)

Tu objetivo: crear tu propio phantom (al menos 2 materiales) y generar su radiografía. Debes incluir 3 figuras: (1) phantom, (2) intensidad del detector, (3) radiografía final ventaneada. Luego escribe un párrafo corto interpretando (5–8 líneas) qué zonas atenúan más y cómo se ven.

Código

def crear_phantom_estudiante(N):
    mu_fondo = 0.02
    phantom = np.full((N, N), mu_fondo, dtype=float)

    # TODO: forma 1 (más atenuante): bloque
    phantom[int(0.30*N):int(0.60*N), int(0.20*N):int(0.35*N)] = 0.09

    # TODO: forma 2 (intermedia): círculo
    y = np.arange(N)[:, None]
    x = np.arange(N)[None, :]
    cy, cx = int(0.65*N), int(0.65*N)
    r = int(0.12*N)
    phantom[(x - cx)**2 + (y - cy)**2 <= r**2] = 0.05

    return phantom

mu_est = crear_phantom_estudiante(N)

# Proyección + detector + radiografía
tau_est = np.sum(mu_est, axis=1) * dx
I_det_est = detector_ideal_polienergetico(tau_est, I0_por_energia)
rad_est, _ = reconstruir_radiografia(I_det_est, I0_total, eps=eps)

# A imagen para visualizar
I_det_est_img = tile_a_imagen(I_det_est, N)
rad_est_img = tile_a_imagen(rad_est, N)

# Ventaneo automático simple
centro_est = float(np.median(rad_est))
ancho_est = float(np.percentile(rad_est, 95) - np.percentile(rad_est, 5))
rad_est_win = ventaneo_lineal(rad_est_img, centro_est, ancho_est)

# 3 figuras requeridas
fig, ax = plt.subplots(1, 3, figsize=(12, 4))

ax[0].imshow(mu_est, cmap="gray")
ax[0].set_title("1) Phantom μ(y,x)")
ax[0].set_axis_off()

ax[1].imshow(I_det_est_img, cmap="gray")
ax[1].set_title("2) Detector ideal I(y) (tiling)")
ax[1].set_axis_off()

ax[2].imshow(rad_est_win, cmap="gray", vmin=0, vmax=1)
ax[2].set_title("3) Radiografía final (ventaneada)")
ax[2].set_axis_off()

plt.tight_layout()
plt.show()

Interpretación (5–8 líneas, escribe aquí en tu entrega): - ¿Dónde está el material más atenuante? ¿Cómo se ve en I_det y en rad? - ¿Qué zonas dejan pasar más fotones? - ¿El ventaneo te ayudó a ver bordes o materiales? ¿Por qué?

12) Entregables (checklist de evaluación)

El código corre sin errores al renderizar el HTML.
Se ven (o se exportan) 3 figuras del mini-proyecto:
- Phantom (mu_est)
- Intensidad en detector (I_det_est)
- Radiografía final (rad_est con ventaneo)
Interpretación corta (5–8 líneas) en lenguaje claro.
(Opcional) Cambia el “espectro simple” (pesos) y comenta qué cambia.
- Pista: con este modelo (sin μ dependiente de energía), rad debería cambiar muy poco.

Apéndice (opcional): exportar figuras a PNG

Si necesitas entregar archivos, usa savefig. No es obligatorio si todo queda visible en el HTML.

Código

# Ejemplo (descomenta si quieres guardar):
# plt.figure()
# plt.imshow(mu_est, cmap="gray")
# plt.axis("off")
# plt.title("Phantom del estudiante")
# plt.savefig("phantom_est.png", dpi=150, bbox_inches="tight")
# plt.close()

--- title: "Laboratorio (hands-on): Reconstrucción básica de una radiografía (Beer–Lambert + haz paralelo)" subtitle: "Fuente → Objeto → Detector → (normalización + log) → Radiografía lista para ver" lang: es author: "Ph.D. Pablo Eduardo Caicedo Rodríguez" date: last-modified format: html: toc: true toc-depth: 3 code-fold: true code-tools: true code-overflow: wrap code-line-numbers: true code-copy: true fig-align: center self-contained: true execute: echo: true warning: false message: false jupyter: python3 --- ## 1) Introducción En este laboratorio vas a construir una “radiografía” sintética desde cero, entendiendo la cadena **fuente → objeto → detector**. La idea central es simple: **si un rayo atraviesa más material (o material más atenuante), llegan menos fotones al detector**. ### ¿Qué entra y qué sale? - **Entra (fuente):** cuántos fotones salen (lo llamaremos `I0`). - **Entra (objeto):** un mapa 2D `mu[y,x]` que dice qué tan atenuante es cada zona. - **Sale (detector):** cuántos fotones llegan por rayo (lo llamaremos `I_det`). - **Sale (radiografía):** una imagen final `rad = -log(I_det/I0)` que es proporcional a “cuánta atenuación acumuló” cada rayo. ### Objetivos de aprendizaje (al finalizar deberías poder) 1) Explicar con palabras simples Beer–Lambert. 2) Construir phantoms 2D con regiones de distintos materiales. 3) Simular proyecciones con **haz paralelo** (sumas en filas). 4) Interpretar `I_det` (más fotones vs menos fotones) y `rad` (más atenuación vs menos). ### Mapa del pipeline (una sola vista) ``` I0 (fuente) → tau (suma de μ en el camino) → I_det = I0·exp(-tau) → rad = -log(I_det/I0) ``` ### Glosario rápido (lo mínimo para no perderse) | Símbolo/variable | Qué es | Cómo se interpreta | |---|---|---| | `mu` (μ) | Mapa 2D del objeto | μ grande = material más atenuante | | `tau` (τ) | Atenuación acumulada por rayo | τ grande = el rayo vio “más bloqueo” | | `I0` | Intensidad sin objeto | referencia (brillo de la fuente) | | `I_det` (I) | Intensidad medida con objeto | baja cuando τ sube | | `rad` | `-log(I_det/I0)` | sube cuando τ sube (radiografía) | ::: {.callout-note} **Regla mental que usarás todo el tiempo** Más `mu` (o más “recorrido” dentro del objeto) → mayor atenuación acumulada `tau` → menor intensidad en el detector `I_det` → mayor valor en la radiografía `rad` (después del log). ::: --- ## 2) Preparación (imports y parámetros) Trabajaremos con arreglos 2D `N×N`. En 2D, el detector “natural” es una señal 1D (un valor por fila). Para visualizarlo como imagen, repetiremos esa señal (tiling) y la mostraremos como 2D. ```{python} import numpy as np import matplotlib.pyplot as plt plt.rcParams["figure.dpi"] = 120 N = 256 # tamaño del phantom (rápido en laptop) dx = 1.0 # tamaño de pixel (unidades arbitrarias) eps = 1e-12 # evita log(0) si hay valores muy pequeños ``` --- ## 3) Modelo físico: Beer–Lambert (muy claro y sin extras) ### 3.1 La idea (en palabras) Imagina el objeto como una **cortina** entre la fuente y el detector. - Cortina “ligera” (μ pequeño) → pasan más fotones. - Cortina “pesada” (μ grande) → pasan menos fotones. `μ(x,y)` es un número que indica qué tan “bloqueador” es el material en ese punto. Un rayo no ve un solo punto: atraviesa una **ruta** y va acumulando bloqueo. ### 3.2 La ecuación $$ I = I_0 \, \exp\left(-\int \mu(x,y)\,dl\right) $$ **Traducción directa (sin misterio):** - `I0`: fotones que llegarían si NO hubiera objeto. - `I`: fotones que llegan DESPUÉS de atravesar el objeto. - `∫ μ dl` (lo llamaremos `tau`): “bloqueo total” acumulado a lo largo del camino del rayo. ### 3.3 Por qué luego usamos log Si definimos $$ \tau = \int \mu\,dl, $$ entonces $$ \frac{I}{I_0} = e^{-\tau} \quad\Rightarrow\quad -\log\left(\frac{I}{I_0}\right)=\tau. $$ O sea: el **log** “deshace” la exponencial y nos devuelve algo proporcional a la atenuación acumulada. Eso es lo que aquí llamamos “radiografía reconstruida”. ::: {.callout-warning} **Qué NO incluye este laboratorio (a propósito):** no hay scatter, no hay ruido Poisson, no hay endurecimiento de haz, no hay blur del detector. Esto es para aprender el esqueleto del pipeline sin distracciones. ::: --- ## 4) Construcción de un phantom (mapa 2D de materiales μ) ### 4.1 Qué es un phantom (didáctico) Un *phantom* es un “objeto inventado” para practicar. Aquí es una imagen 2D de números `mu[y,x]`: cada pixel dice “qué material hay” (qué tanto atenúa). - `mu` pequeño → zona más “transparente” (deja pasar más fotones). - `mu` grande → zona más “opaca” (deja pasar menos fotones). El fondo representa un material base (`mu_fondo`). Las formas internas tienen otro μ para crear **contraste**. ### 4.2 Cómo “dibujamos” materiales con máscaras En NumPy podemos construir una **máscara booleana** (True/False). - True = “píntalo con el material del círculo/bloque”. - False = “déjalo como fondo”. En el siguiente bloque creamos dos phantoms: A) fondo + círculo, B) fondo + bloque + círculo. ```{python} def crear_phantom_circulo(N, mu_fondo=0.02, mu_circulo=0.06, radio_frac=0.25, centro=None): """Phantom A: fondo uniforme + un círculo.""" if centro is None: centro = (N//2, N//2) # (y, x) # Grid de coordenadas y = np.arange(N)[:, None] x = np.arange(N)[None, :] cy, cx = centro r = radio_frac * N # 1) Empieza con fondo phantom = np.full((N, N), mu_fondo, dtype=float) # 2) Máscara del círculo (True dentro del círculo) mascara = (x - cx)**2 + (y - cy)**2 <= r**2 # 3) "Pinta" el círculo cambiando mu phantom[mascara] = mu_circulo return phantom def crear_phantom_dos_materiales(N, mu_fondo=0.02, mu_bloque=0.10, mu_circulo=0.05): """Phantom B: fondo + bloque (más atenuante) + círculo (intermedio).""" phantom = np.full((N, N), mu_fondo, dtype=float) # Bloque rectangular (material más atenuante) y1, y2 = int(0.35*N), int(0.75*N) x1, x2 = int(0.15*N), int(0.45*N) phantom[y1:y2, x1:x2] = mu_bloque # Círculo (material intermedio) y = np.arange(N)[:, None] x = np.arange(N)[None, :] cy, cx = int(0.45*N), int(0.70*N) r = int(0.15*N) mascara = (x - cx)**2 + (y - cy)**2 <= r**2 phantom[mascara] = mu_circulo return phantom mu_A = crear_phantom_circulo(N) mu_B = crear_phantom_dos_materiales(N) fig, ax = plt.subplots(1, 2, figsize=(9, 4)) im0 = ax[0].imshow(mu_A, cmap="gray") ax[0].set_title("Phantom A: círculo en fondo (mapa μ)") plt.colorbar(im0, ax=ax[0], fraction=0.046) im1 = ax[1].imshow(mu_B, cmap="gray") ax[1].set_title("Phantom B: dos materiales (mapa μ)") plt.colorbar(im1, ax=ax[1], fraction=0.046) for a in ax: a.set_axis_off() plt.tight_layout() plt.show() ``` ::: {.callout-tip} **Lectura rápida del phantom** En la imagen del phantom NO estás viendo una radiografía. Estás viendo “materiales”: dónde μ es mayor, ese material atenúa más. ::: ### Ejercicio guiado 1 (predicción antes de correr el pipeline) 1) Si aumentas `mu_circulo`, ¿la intensidad `I_det` en esas filas sube o baja? ¿por qué? 2) Si aumentas `mu_bloque`, ¿qué esperas que pase con `I_det` y con `rad`? Escribe tu predicción en 2–3 líneas. Luego la comprobarás cuando veas `I_det` y `rad`. --- ## 5) Generador polienergético idealizado (espectro simple) ### 5.1 Qué significa “polienergético” aquí En rayos X reales, la fuente emite fotones con diferentes energías (un espectro). Aquí lo simplificamos a **3 energías discretas** con **pesos**. Esto es didáctico. **Importante:** en este laboratorio NO hacemos que μ dependa de energía. Por eso, el “espectro” se reduce a un `I0_total` efectivo, pero nos sirve para entender la idea de “sumar contribuciones”. ### 5.2 Implementación (3 energías + pesos) ```{python} energias_keV = np.array([30, 60, 90], dtype=float) pesos = np.array([0.50, 0.30, 0.20], dtype=float) pesos = pesos / pesos.sum() I0_total = 1e5 I0_por_energia = I0_total * pesos print("Energías (keV):", energias_keV) print("Pesos:", pesos) print("I0 por energía:", I0_por_energia) print("I0 efectivo:", I0_por_energia.sum()) # Visualización del espectro discreto (solo para ver pesos) fig, ax = plt.subplots(figsize=(5.5, 3.2)) ax.bar(energias_keV.astype(int), pesos) ax.set_title("Espectro simple (discreto): pesos por energía") ax.set_xlabel("Energía (keV)") ax.set_ylabel("Peso (suma = 1)") plt.tight_layout() plt.show() ``` **Pregunta corta (1 línea):** Si μ no cambia con energía, ¿por qué esperarías que `rad = -log(I/I0)` cambie muy poco al cambiar `pesos`? --- ## 6) Simulación de proyecciones (haz paralelo, integral discreta) ### 6.1 Qué es “proyectar” en este laboratorio Con **haz paralelo** y en 2D, usaremos rayos horizontales. Para cada fila `y`, el rayo atraviesa todos los pixeles `x` y acumula atenuación: $$ \tau(y) \approx \sum_x \mu(y,x)\,\Delta x $$ `tau(y)` es una señal 1D: un valor por fila (un valor por rayo). ### 6.2 Implementación (sumar en x) ```{python} def proyectar_haz_paralelo(mu, dx=1.0): """Proyección 2D→1D: tau[y] = sum_x mu[y,x]*dx""" tau = np.sum(mu, axis=1) * dx return tau tau_A = proyectar_haz_paralelo(mu_A, dx=dx) tau_B = proyectar_haz_paralelo(mu_B, dx=dx) fig, ax = plt.subplots(1, 2, figsize=(9, 3.2)) ax[0].plot(tau_A) ax[0].set_title("Camino óptico τ(y) - A") ax[0].set_xlabel("índice y (pixel)") ax[0].set_ylabel("τ (aprox.)") ax[1].plot(tau_B) ax[1].set_title("Camino óptico τ(y) - B") ax[1].set_xlabel("índice y (pixel)") plt.tight_layout() plt.show() ``` ### 6.2b Conectar “una fila del phantom” con un valor de τ Para que no se vuelva abstracto: el rayo de la fila `y0` atraviesa todos los pixeles de esa fila. Entonces `tau[y0]` es literalmente la suma de `mu[y0, :]` (por `dx`). En el siguiente bloque miramos una fila y verificamos la suma. ```{python} # Escoge una fila (rayo) para inspeccionar y0 = N // 2 # Perfil de mu en esa fila perfil_A = mu_A[y0, :] perfil_B = mu_B[y0, :] print(f"Fila y0={y0}: tau_A[y0]={tau_A[y0]:.4f} (suma discreta)") print(f"Fila y0={y0}: tau_B[y0]={tau_B[y0]:.4f} (suma discreta)") fig, ax = plt.subplots(1, 2, figsize=(9, 3)) ax[0].plot(perfil_A) ax[0].set_title("Perfil μ en una fila (A)") ax[0].set_xlabel("x (pixel)") ax[0].set_ylabel("μ") ax[1].plot(perfil_B) ax[1].set_title("Perfil μ en una fila (B)") ax[1].set_xlabel("x (pixel)") plt.tight_layout() plt.show() ``` ::: {.callout-note} **Interpretación de τ(y)** Si en cierta fila el rayo cruza el círculo o el bloque, τ(y) sube (porque suma μ más grande). Más τ → menos intensidad en el detector (por la exponencial). ::: ### 6.3 Para ver 1D como “imagen” (solo visual) Esto NO agrega física: solo repite la señal 1D en columnas para verla como imagen. ```{python} def tile_a_imagen(signal_1d, N): """Convierte una señal 1D (N,) a una imagen (N,N) repitiéndola en x.""" return np.tile(signal_1d[:, None], (1, N)) tau_A_img = tile_a_imagen(tau_A, N) tau_B_img = tile_a_imagen(tau_B, N) fig, ax = plt.subplots(1, 2, figsize=(9, 4)) im0 = ax[0].imshow(tau_A_img, cmap="magma") ax[0].set_title("τ(y) como imagen (tiling) - A") plt.colorbar(im0, ax=ax[0], fraction=0.046) im1 = ax[1].imshow(tau_B_img, cmap="magma") ax[1].set_title("τ(y) como imagen (tiling) - B") plt.colorbar(im1, ax=ax[1], fraction=0.046) for a in ax: a.set_axis_off() plt.tight_layout() plt.show() ``` --- ## 7) Detector ideal (sin blur, sin ruido) ### 7.1 Qué calcula el detector El detector ideal “cuenta” fotones que llegan por rayo. Para cada energía `k`: $$ I_k(y) = I0_k\,e^{-\tau(y)} $$ y luego sumamos todas las energías: $$ I(y) = \sum_k I_k(y). $$ ### 7.2 Implementación ```{python} def detector_ideal_polienergetico(tau, I0_por_energia): """I(y) = sum_k I0_k * exp(-tau(y))""" I = np.zeros_like(tau, dtype=float) for I0_k in I0_por_energia: I += I0_k * np.exp(-tau) return I I_det_A = detector_ideal_polienergetico(tau_A, I0_por_energia) I_det_B = detector_ideal_polienergetico(tau_B, I0_por_energia) fig, ax = plt.subplots(1, 2, figsize=(9, 3.2)) ax[0].plot(I_det_A) ax[0].set_title("Intensidad ideal en detector I(y) - A") ax[0].set_xlabel("índice y (pixel)") ax[0].set_ylabel("conteos (ideal)") ax[1].plot(I_det_B) ax[1].set_title("Intensidad ideal en detector I(y) - B") ax[1].set_xlabel("índice y (pixel)") plt.tight_layout() plt.show() ``` ### 7.2b Ver la exponencial con tus ojos (τ vs I) Si Beer–Lambert es correcto, al aumentar `tau` la intensidad debe caer de forma exponencial. En este bloque graficamos `I_det` contra `tau` para ver esa caída. ```{python} fig, ax = plt.subplots(figsize=(5.8, 4)) ax.scatter(tau_A, I_det_A, s=10, alpha=0.6, label="A") ax.scatter(tau_B, I_det_B, s=10, alpha=0.6, label="B") ax.set_title("Relación exponencial: I_det vs τ") ax.set_xlabel("τ (atenuación acumulada)") ax.set_ylabel("I_det (conteos ideales)") ax.legend() plt.tight_layout() plt.show() ``` ::: {.callout-tip} **Cómo “leer” I_det** Picos altos = llegaron muchos fotones (poca atenuación). Valles bajos = llegaron pocos fotones (mucha atenuación). ::: ### Ejercicio guiado 2 (comprobación de tu predicción) Vuelve al Ejercicio 1. Ahora que ves `I_det`, confirma si tu predicción fue correcta y escribe 2 líneas: - “Subí μ en X → I_det en esa zona **bajó/subió** porque …”. --- ## 8) Reconstrucción tipo radiografía: normalización + log ### 8.1 Por qué hacemos estos dos pasos 1) **Normalización:** `I_det/I0` elimina el “brillo” de la fuente y deja solo el efecto del objeto. 2) **Log negativo:** como `I/I0 = exp(-tau)`, entonces `-log(I/I0) = tau`. Esto es clave: la radiografía reconstruida `rad` es proporcional al “camino óptico” `tau`. ### 8.2 Implementación ```{python} def reconstruir_radiografia(I_det, I0_total, eps=1e-12): """rad(y) = -log( I_det(y) / I0_total )""" I_norm = (I_det + eps) / (I0_total + eps) rad = -np.log(I_norm) return rad, I_norm rad_A, I_norm_A = reconstruir_radiografia(I_det_A, I0_total, eps=eps) rad_B, I_norm_B = reconstruir_radiografia(I_det_B, I0_total, eps=eps) # Verificación visual: rad ~ tau fig, ax = plt.subplots(1, 2, figsize=(9, 3.2)) ax[0].plot(tau_A, label="tau (verdad)") ax[0].plot(rad_A, "--", label="rad (-log)") ax[0].set_title("A: tau vs rad") ax[0].legend() ax[1].plot(tau_B, label="tau (verdad)") ax[1].plot(rad_B, "--", label="rad (-log)") ax[1].set_title("B: tau vs rad") ax[1].legend() plt.tight_layout() plt.show() ``` ::: {.callout-note} **Interpretación de rad** Valores altos en `rad` significan: el rayo vio más material (o material más atenuante). Esto suele ser lo que queremos “ver” en una radiografía. ::: ### Ejercicio guiado 3 (invarianza a I0) Multiplica `I0_total` por 10 (por ejemplo `1e6`) y repite el pipeline. - ¿`I_det` cambia? (sí) - ¿`rad` cambia mucho? (debería cambiar muy poco) --- ## 9) Visualización: ventaneo (contraste simple) ### 9.1 Por qué ventaneamos Una radiografía es una matriz de valores, pero el ojo necesita un rango limitado para ver detalles. Ventaneo = escoger un intervalo `[low, high]` y mapearlo a [0,1]. - Ancho pequeño → más contraste, pero se saturan zonas. - Ancho grande → menos contraste, más tonos visibles. ### 9.2 Implementación ```{python} def ventaneo_lineal(img, centro, ancho): """ Ventaneo simple: - recorta a [centro - ancho/2, centro + ancho/2] - reescala a [0, 1] """ low = centro - ancho/2 high = centro + ancho/2 out = np.clip(img, low, high) out = (out - low) / (high - low + 1e-12) return out # Convertimos señales 1D a imágenes 2D (tiling) para visualizar rad_A_img = tile_a_imagen(rad_A, N) rad_B_img = tile_a_imagen(rad_B, N) # Ventaneo automático simple con percentiles centro_A = float(np.median(rad_A)) ancho_A = float(np.percentile(rad_A, 95) - np.percentile(rad_A, 5)) centro_B = float(np.median(rad_B)) ancho_B = float(np.percentile(rad_B, 95) - np.percentile(rad_B, 5)) rad_A_win = ventaneo_lineal(rad_A_img, centro_A, ancho_A) rad_B_win = ventaneo_lineal(rad_B_img, centro_B, ancho_B) fig, ax = plt.subplots(1, 2, figsize=(9, 4)) ax[0].imshow(rad_A_win, cmap="gray", vmin=0, vmax=1) ax[0].set_title("Radiografía A (ventaneada)") ax[1].imshow(rad_B_win, cmap="gray", vmin=0, vmax=1) ax[1].set_title("Radiografía B (ventaneada)") for a in ax: a.set_axis_off() plt.tight_layout() plt.show() ``` ### 9.3 ¿Qué está haciendo el ventaneo? (histograma) El histograma te muestra cuántos pixeles tienen cierto valor de `rad`. Al ventanear, elegimos un rango `[low, high]` y todo lo demás se recorta. Aquí dibujamos ese rango sobre el histograma para entenderlo. ```{python} low_A = centro_A - ancho_A/2 high_A = centro_A + ancho_A/2 fig, ax = plt.subplots(figsize=(6.2, 3.6)) ax.hist(rad_A, bins=40) ax.axvline(low_A, linestyle='--', label='low (A)') ax.axvline(high_A, linestyle='--', label='high (A)') ax.set_title("Histograma de rad (A) y rango de ventaneo") ax.set_xlabel("rad") ax.set_ylabel("número de filas (pixeles en 1D)") ax.legend() plt.tight_layout() plt.show() ``` ### Ejercicio guiado 4 (control de contraste) Prueba estos cambios y describe en 2–3 líneas qué notas: - Reduce `ancho_A` a la mitad. ¿Qué detalle aparece y qué se satura? - Aumenta `ancho_A` al doble. ¿Qué se ve más “plano” y por qué? --- ## 10) Comparación final: phantom vs detector vs radiografía Aquí juntamos todo para leer el pipeline de un vistazo: - `mu` = “qué hay” (materiales). - `I_det` = “qué midió el detector” (conteos). - `rad` = “radiografía” (atenuación acumulada). ```{python} def resumen_pipeline(mu, I0_total, I0_por_energia, dx=1.0, eps=1e-12): tau = np.sum(mu, axis=1) * dx I_det = detector_ideal_polienergetico(tau, I0_por_energia) rad, _ = reconstruir_radiografia(I_det, I0_total, eps=eps) return tau, I_det, rad tauA, IA, radA = resumen_pipeline(mu_A, I0_total, I0_por_energia, dx=dx, eps=eps) tauB, IB, radB = resumen_pipeline(mu_B, I0_total, I0_por_energia, dx=dx, eps=eps) fig, ax = plt.subplots(2, 3, figsize=(10, 6)) ax[0,0].imshow(mu_A, cmap="gray"); ax[0,0].set_title("A: phantom μ(y,x)") ax[1,0].imshow(mu_B, cmap="gray"); ax[1,0].set_title("B: phantom μ(y,x)") ax[0,1].imshow(tile_a_imagen(IA, N), cmap="gray"); ax[0,1].set_title("A: detector I(y) (tiling)") ax[1,1].imshow(tile_a_imagen(IB, N), cmap="gray"); ax[1,1].set_title("B: detector I(y) (tiling)") ax[0,2].imshow(tile_a_imagen(radA, N), cmap="inferno"); ax[0,2].set_title("A: rad = -log(I/I0)") ax[1,2].imshow(tile_a_imagen(radB, N), cmap="inferno"); ax[1,2].set_title("B: rad = -log(I/I0)") for a in ax.ravel(): a.set_axis_off() plt.tight_layout() plt.show() ``` --- ## 11) Mini-proyecto final (entregable) Tu objetivo: crear tu propio phantom (al menos 2 materiales) y generar su radiografía. Debes incluir 3 figuras: (1) phantom, (2) intensidad del detector, (3) radiografía final ventaneada. Luego escribe un párrafo corto interpretando (5–8 líneas) qué zonas atenúan más y cómo se ven. ```{python} def crear_phantom_estudiante(N): mu_fondo = 0.02 phantom = np.full((N, N), mu_fondo, dtype=float) # TODO: forma 1 (más atenuante): bloque phantom[int(0.30*N):int(0.60*N), int(0.20*N):int(0.35*N)] = 0.09 # TODO: forma 2 (intermedia): círculo y = np.arange(N)[:, None] x = np.arange(N)[None, :] cy, cx = int(0.65*N), int(0.65*N) r = int(0.12*N) phantom[(x - cx)**2 + (y - cy)**2 <= r**2] = 0.05 return phantom mu_est = crear_phantom_estudiante(N) # Proyección + detector + radiografía tau_est = np.sum(mu_est, axis=1) * dx I_det_est = detector_ideal_polienergetico(tau_est, I0_por_energia) rad_est, _ = reconstruir_radiografia(I_det_est, I0_total, eps=eps) # A imagen para visualizar I_det_est_img = tile_a_imagen(I_det_est, N) rad_est_img = tile_a_imagen(rad_est, N) # Ventaneo automático simple centro_est = float(np.median(rad_est)) ancho_est = float(np.percentile(rad_est, 95) - np.percentile(rad_est, 5)) rad_est_win = ventaneo_lineal(rad_est_img, centro_est, ancho_est) # 3 figuras requeridas fig, ax = plt.subplots(1, 3, figsize=(12, 4)) ax[0].imshow(mu_est, cmap="gray") ax[0].set_title("1) Phantom μ(y,x)") ax[0].set_axis_off() ax[1].imshow(I_det_est_img, cmap="gray") ax[1].set_title("2) Detector ideal I(y) (tiling)") ax[1].set_axis_off() ax[2].imshow(rad_est_win, cmap="gray", vmin=0, vmax=1) ax[2].set_title("3) Radiografía final (ventaneada)") ax[2].set_axis_off() plt.tight_layout() plt.show() ``` **Interpretación (5–8 líneas, escribe aquí en tu entrega):** - ¿Dónde está el material más atenuante? ¿Cómo se ve en `I_det` y en `rad`? - ¿Qué zonas dejan pasar más fotones? - ¿El ventaneo te ayudó a ver bordes o materiales? ¿Por qué? --- ## 12) Entregables (checklist de evaluación) - [ ] El código corre sin errores al renderizar el HTML. - [ ] Se ven (o se exportan) 3 figuras del mini-proyecto: - [ ] Phantom (`mu_est`) - [ ] Intensidad en detector (`I_det_est`) - [ ] Radiografía final (`rad_est` con ventaneo) - [ ] Interpretación corta (5–8 líneas) en lenguaje claro. - [ ] (Opcional) Cambia el “espectro simple” (`pesos`) y comenta qué cambia. - Pista: con este modelo (sin μ dependiente de energía), `rad` debería cambiar muy poco. --- ## Apéndice (opcional): exportar figuras a PNG Si necesitas entregar archivos, usa `savefig`. No es obligatorio si todo queda visible en el HTML. ```{python} # Ejemplo (descomenta si quieres guardar): # plt.figure() # plt.imshow(mu_est, cmap="gray") # plt.axis("off") # plt.title("Phantom del estudiante") # plt.savefig("phantom_est.png", dpi=150, bbox_inches="tight") # plt.close() ```